function [ResEuler,yGrid,xGrid] = EulerEqError_demo(model,extendedPerOn,orderApp,setup)
% Some indices
ny      = model.ny;        %Number of output variables
nx      = model.nx;        %Number of state variables

%Unfolding parameters
names = fieldnames(model.params);
for i=1:length(names)
    eval([names{i},'=',num2str(model.params.(names{i}),16),';']);
end

% The raw grid for the states
if setup.xPointsDim == 0
    xGrid = setup.xSim';
else
    dx       = (setup.upperx-setup.lowerx)/(setup.xPointsDim-1);
    xGridtmp = zeros(nx,setup.xPointsDim);
    for j=1:setup.xPointsDim
        xGridtmp(:,j) = setup.lowerx(:,1)+(j-1)*dx(:,1);
    end
    tmp = cell(1,nx);
    for i=1:nx
        tmp(i) = mat2cell(xGridtmp(i,:),1,setup.xPointsDim);
    end
    xGrid = cartprod(tmp);
end

% The grid for numerical itegration
tmp_e        = cell(1,model.ne);
tmp_w        = cell(1,model.ne);
for i=1:model.ne
    tmp_e(i) = mat2cell(setup.GH_e,length(setup.GH_e),1);
    tmp_w(i) = mat2cell(setup.GH_w,length(setup.GH_w),1);
end
eGrid = cartprod(tmp_e);
wGrid = prod(cartprod(tmp_w),2);
%% Building the residuals
f        = zeros(ny+nx,1);
ResEuler = zeros(size(xGrid,1),ny+nx);
yGrid    = zeros(size(xGrid,1),ny);
version  = 1*0;
for k = 1:size(xGrid,1)
    %k
    x_cu   = xGrid(k,:)';
    % Getting the states at time t
    k_cu  = model.h0(1,1) + x_cu(1,1);
    c_ba1 = model.h0(2,1) + x_cu(2,1);
    a_cu  = model.h0(3,1) + x_cu(3,1);
    d_cu  = model.h0(4,1) + x_cu(4,1);
    
    % Getting the controls at time t and part of the states at time t+1
    if extendedPerOn == 1
        [yVec_cu,xVec_cup,diagOut] = policyFunctionEPer(x_cu+model.h0,model,setup.setupEPer);
    else
        [yVec_cu,xVec_cup]=policy(x_cu+model.h0,model.g0,model.h0,model.derivs,orderApp);
    end
    yGrid(k,:) = (yVec_cu-model.g0)';
    % The controls
    c_cu  = yVec_cu(1,1);
    iv_cu = yVec_cu(2,1);
    y_cu  = yVec_cu(3,1);
    la_cu = yVec_cu(4,1);
    n_cu  = yVec_cu(5,1);
    rk_cu = yVec_cu(6,1);
    w_cu  = yVec_cu(7,1);
    
    Euler = zeros(ny+nx,1);
    if version == 1
        for i=1:length(wGrid)
            x_cup = xVec_cup-model.h0;
            % Adding shocks to the states
            x_cup = x_cup + model.eta*eGrid(i,:)';
            
            % Getting controls at time t+1
            if extendedPerOn == 1
                yVec_cup = policyFunctionEPer(x_cup+model.h0,model,setup.setupEPer);
            else
                yVec_cup = policy(x_cup+model.h0,model.g0,model.h0,model.derivs,orderApp);
            end
            c_cup  = yVec_cup(1,1);
            iv_cup = yVec_cup(2,1);
            y_cup  = yVec_cup(3,1);
            la_cup = yVec_cup(4,1);
            n_cup  = yVec_cup(5,1);
            rk_cup = yVec_cup(6,1);
            w_cup  = yVec_cup(7,1);
            
            % Getting the states at time t+1
            %for i=1:nx
            %    eval([char(model.xp(i)),'=',num2str(x_cup(i,1),16),';']);
            %end
            k_cup  = model.h0(1,1) + x_cup(1,1);
            c_ba1p = model.h0(2,1) + x_cup(2,1);
            a_cup  = model.h0(3,1) + x_cup(3,1);
            d_cup  = model.h0(4,1) + x_cup(4,1);
            
            %% Model Specific
            % START DISPLAYING f
            f(1,1)=exp(-la_cu)*exp(d_cu)*(1/(exp(c_cu) - B*exp(c_ba1))^ETAc - (B*BETTA*exp(d_cup))/(exp(c_cup) - B*exp(c_cu))^ETAc) - 1;
            f(2,1)=1 - (THETA*exp(-la_cu)*exp(-w_cu)*exp(d_cu))/(1 - exp(n_cu))^ETAl;
            f(3,1)=BETTA*exp(-la_cu)*exp(la_cup)*(exp(rk_cup) - DELTA + 1) - 1;
            f(4,1)=exp(rk_cu) + (exp(a_cu)*exp(n_cu)^ALFA*(ALFA - 1))/exp(k_cu)^ALFA;
            f(5,1)=1 - ALFA*exp(-w_cu)*exp(a_cu)*exp(k_cu)^(1 - ALFA)*exp(n_cu)^(ALFA - 1);
            f(6,1)=1 - exp(-y_cu)*(exp(c_cu) + exp(iv_cu));
            f(7,1)=exp(-y_cu)*exp(a_cu)*exp(k_cu)^(1 - ALFA)*exp(n_cu)^ALFA - 1;
            f(8,1)=1 - exp(-iv_cu)*(exp(k_cup) + exp(k_cu)*(DELTA - 1));
            f(9,1)=1 - exp(-c_cu)*exp(c_ba1p);
            f(10,1)=RHOA*log(exp(a_cu)) - log(exp(a_cup));
            f(11,1)=RHOD*log(exp(d_cu)) - log(exp(d_cup));
            % END DISPLAYING f
            Euler(:,1) = Euler(:,1) + wGrid(i,1)*f(:,1);
        end
    else
        for i1_eps=1:size(setup.GH_e,1)
            for i2_eps=1:size(setup.GH_e,1)
                %% Model Specific
                x_cup= xVec_cup-model.h0;
                % Adding shocks to the states
                % For a_cu
                x_cup(3,1) = x_cup(3,1) + STDA*setup.GH_e(i1_eps,1);
                % for d_cu
                x_cup(4,1) = x_cup(4,1) + STDD*setup.GH_e(i2_eps,1);
                
                %% Getting controls at time t+1
                if extendedPerOn == 1
                    yVec_cup = policyFunctionEPer(x_cup+model.h0,model,setup.setupEPer);
                else
                    yVec_cup =policy(x_cup+model.h0,model.g0,model.h0,model.derivs,orderApp);
                end
                c_cup  = yVec_cup(1,1);
                iv_cup = yVec_cup(2,1);
                y_cup  = yVec_cup(3,1);
                la_cup = yVec_cup(4,1);
                n_cup  = yVec_cup(5,1);
                rk_cup = yVec_cup(6,1);
                w_cup  = yVec_cup(7,1);
                
                % Getting the states at time t+1
                %for i=1:nx
                %    eval([char(model.xp(i)),'=',num2str(x_cup(i,1),16),';']);
                %end
                k_cup  = model.h0(1,1) + x_cup(1,1);
                c_ba1p = model.h0(2,1) + x_cup(2,1);
                a_cup  = model.h0(3,1) + x_cup(3,1);
                d_cup  = model.h0(4,1) + x_cup(4,1);
                
                %% Model Specific
                % START DISPLAYING f
                f(1,1)=exp(-la_cu)*exp(d_cu)*(1/(exp(c_cu) - B*exp(c_ba1))^ETAc - (B*BETTA*exp(d_cup))/(exp(c_cup) - B*exp(c_cu))^ETAc) - 1;
                f(2,1)=1 - (THETA*exp(-la_cu)*exp(-w_cu)*exp(d_cu))/(1 - exp(n_cu))^ETAl;
                f(3,1)=BETTA*exp(-la_cu)*exp(la_cup)*(exp(rk_cup) - DELTA + 1) - 1;
                f(4,1)=exp(rk_cu) + (exp(a_cu)*exp(n_cu)^ALFA*(ALFA - 1))/exp(k_cu)^ALFA;
                f(5,1)=1 - ALFA*exp(-w_cu)*exp(a_cu)*exp(k_cu)^(1 - ALFA)*exp(n_cu)^(ALFA - 1);
                f(6,1)=1 - exp(-y_cu)*(exp(c_cu) + exp(iv_cu));
                f(7,1)=exp(-y_cu)*exp(a_cu)*exp(k_cu)^(1 - ALFA)*exp(n_cu)^ALFA - 1;
                f(8,1)=1 - exp(-iv_cu)*(exp(k_cup) + exp(k_cu)*(DELTA - 1));
                f(9,1)=1 - exp(-c_cu)*exp(c_ba1p);
                f(10,1)=RHOA*log(exp(a_cu)) - log(exp(a_cup));
                f(11,1)=RHOD*log(exp(d_cu)) - log(exp(d_cup));
                % END DISPLAYING f
                Euler(:,1) = Euler(:,1) + setup.GH_w(i2_eps)*setup.GH_w(i1_eps)*f(:,1);
            end
        end
    end
    ResEuler(k,:) = Euler(:,1)';
end

end